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Abstract. This work develops novel error expansions with computable lead- 
ing order terms for the global weak error in the tau-leap discretization of pure 
jump processes arising in kinetic Monte Carlo models. Accurate computable 
a posteriori error approximations are the basis for adaptive algorithms; a fun- 
damental tool for numerical simulation of both deterministic and stochastic 
dynamical systems. These pure jump processes are simulated either by the 
tau-leap method, or by exact simulation, also referred to as dynamic Monte 
Carlo, the Gillespie algorithm or the Stochastic simulation algorithm. Two 
types of estimates are presented: an a priori estimate for the relative error 
that gives a comparison between the work for the two methods depending on 
the propensity regime, and an a posteriori estimate with computable leading 
order term. 



In this work we derive a global weak error expansion with computable leading 
order term for the tau-leap method. The tau-leap method was originally proposed 
by Gillespie in [10], for approximating homogeneous and well stirred stochastically 
reacting chemical systems. In such systems different species undergo reactions at 
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random times to form new species or to decay. Each reaction can be modeled by 
a propensity junction, directly related to the number of particles, indicating the 
probability of a reaction to happen per unit of time. The notion well stirred here 
means that the number of reactive particle collisions are low compared to the total 
number of collisions. Depending on the number of particles in the system and the 
time to next reaction the reaction process can be modeled differently: 

For all possible regimes, a homogeneous well stirred system in thermal equi- 
librium can be modeled by the chemical master equation (CME). This ordinary 
differential equation describes the time evolution of the probability of each particle 
configuration, cf. [ ] . Since the dimension of the solution space is of the order of 
all the possible configurations, the CME is in practice impossible to solve. On the 
other hand, the system can still be simulated exactly using the Stochastic simulation 
algorithm (SSA), also introduced by Gillespie in [f 1]. The SSA numerically sim- 
ulates the Markov process described by the CME by using dynamic Monte Carlo 
sampling, also referred to as kinetic Monte Carlo. Although the SSA generates 
exact path realizations for the Markov process it is only tractable for low propensi- 
ties. Indeed, since each reaction is simulated exactly by sampling the next reaction 
to happen and the time to this reaction, the total computational work becomes 
roughly inversely proportional to the total propensity. The tau-leap method, on 
the other hand, approximates the SSA by evolving the chemical system with fixed 
time steps, keeping the propensity fixed in each time step, and can be seen as a for- 
ward Euler method for a stochastic differential equation driven by Poisson random 
measures, cf. [ ]. In the limit, as the time steps go to zero, the tau-leap solution 
converges to the SSA, see [24]. 

As the number of particles in the system grows the SSA is sometimes approxi- 
mated by the chemical Langevin diffusion equation. Further, if both the number of 
particles and the total volume of the system goes to infinity at the same speed, i.e. 
the randomness of the system becomes negligible, then the concentration of each 
species over time can be modeled by a deterministic system of ordinary differen- 
tial equations, the reaction rate equations (RREs). Because of different time scales 
in the reactions, RREs are often are stiff, indicating the need for stable tau-leap 
methods and adaptive time-stepping. Implicit stable tau-leap methods that deal 
with stiff cases have been proposed by [4, 23]. Several authors have discussed the 
importance of leap selection procedures to increase efficiency [6, 12] and to avoid 
negative populations [f, 5, 7, 26]. 

Representing the number of particles for the species in the system by the stochas- 
tic vector X(t), the goal in this work is to approximate the real valued quantity, 
E[g(A(T))], for some given function g and initial configuration X(0). We here 
derive an a posteriori estimate for the global weak error E[g(A(T)) — g{X{T))\ 
between exact (SSA) solution X and the tau-leap solution X, based on an exact 
global error representation, cf. Lemma 4.1. The error representation uses the value 
function defined by the Kolmogorov backward equation (2.1). Also, an a priori 
estimate for the relative global weak error that is independent of the number of 
particles is derived. Both the a priori and a posteriori error estimates are based on 
a continuous extension and point wise bounds of the discretely defined value func- 
tion and its derivatives. In [16], similar L-infinity bounds, exponentially dependent 
on the propensity function, were derived. This work improves those estimates show- 
ing that the value function and its derivatives only grows polynomially with the 
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population. This result is based on weighted estimates and a stochastic represen- 
tation of the weighted derivatives in terms of pure jump processes with modified 
propensities. 

Adaptive time stepping algorithms based on computable a posteriori error es- 
timates are fundamental tools for numerical simulation of both deterministic and 
stochastic dynamical systems. In our case, the leading order term of the error ex- 
pansion is approximated by a discrete dual weighted propensity residual, similar 
to [18, 19] for deterministic differential equations, and [20, 21, 22, 25] for stochas- 
tic differential equations of diffusion and jump-diffusion type. To the best of the 
authors' knowledge there are still no results on optimal adaptive time stepping 
algorithms based on a posteriori error estimates for the global weak error in the 
tau-leap method. 

The outline of this work is as follows: first, Section 2 presents the mathematical 
setting together with the main assumptions and the relation with the Kolmogorov 
backward equation upon which our results are based. Next, Section 3 introduces 
the tau-leap method and a numerical scheme for avoiding negative populations 
(following [..]), here called the Poisson bridge tau-leap method. Section 4 presents 
the main results, namely an a priori bound on the relative error of the tau-leap 
method (Theorem 4.4), and a computable a posteriori error approximation using 
a discrete dual (Theorem 4.15). Also, from the a priori bound in Theorem 4.4 a 
work comparison between the tau-leap method and SSA is made. Finally, Section 
5 provides a numerical verification of our a posteriori error estimate in Theorem 
4.15. 

2. Problem 

We consider a well stirred system of d chemical species, which interact through 
M chemical reaction channels. The system is assumed to be confined to a constant 
volume £1 and to be in thermal, but not necessarily chemical, equilibrium at some 
constant temperature. With JQ denoting the number of molecules of species 
i in the system at time t, we want to study the evolution of the state vector 
X t — {X[ ; , . . . , X^), given that the system was initially in some state X to — xq. 
It is here assumed that X t ,xo £ where Z + denotes the set of non- negative 
integers. The goal of the computation is the approximation of the quantity 

E[.9(A T )], 

where g : R d — > R is a given function. 

Each of the j — 1, . . . , M reactions is characterized by the change 

X — > X + Vj 

where Vj £ Z d is a stoichiometric vector that indicates the change in the state 
vector x € Zi, produced by a single firing of the reaction j. Now that we have 
defined what happens when each reaction takes place, we also need to indicate how 
often each of those reactions may occur in time. This information is carried by the 
propensity functions, a,j : — > R, j = 1, . . . , M, which tell the probability of the 
reaction j happening during the infinitesimal interval (t,t 4- dt) i.e. 

P ( j fires during (t,t + dt) | Xt = x j = a,j (x)dt, j = 1, . . . , M. 
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Usually the functions a.,- are just polynomials, see Remark 2.2. The compensator 
of our reaction process is state dependent, and the Kolmogorov backward equation 
for this pure jump process and a smooth observable g : R d -> K is 

M 

t)-u(jc,i)) =0, inZ^_ x [0,T), 

(2-1) j=i 

u = g, on %\ x {T}, 

with the corresponding value function 

u(x,t) :=E[g(X T ) \X t = x\. 

Here, natural boundary conditions restrict the species to the set Zl, see Remark 
2.3. 

Let Z := {i/i, . . . , I'm}- It is interesting to see that the previous equation can be 
rewritten in the following way 

d t u + a (x) (u(x + z,t) — u(x,t) )q(x,dz) = 0, inZlx[0,T), 
(2.2) Jz y ' 

u(-,T)=g, onZ^x{T}, 

M 

a o( x ) ■= J2 a ^ x ^ 

i=i 

where the compensator measure q(x, dz) is atomistic in the stoichiometric vector 
set Z and such that 

q(x, z = v 5 ) = P(z - Vj \X t = x) = j = 1, . . . , M. 

a Q (x) 

Although (2.1) and (2.2) are equivalent they may lead to different discretizations. 
In particular, the stochastic simulation algorithm (SSA) introduced by Gillespie 
[11] is the natural discretization scheme for (2.2). Moreover, the corresponding 
error estimates and error control algorithms behave differently. 

Since (2.1) is driven by pure jumps with finite activity it is in principle possible 
to simulate trajectories of X t exactly; the SSA does precisely that. It only requires 
the sample of two random variables per time step: one to find the time of the next 
reaction and another to decide which is the reaction that is firing at that time. 
The drawback of this algorithm appears clearly as the sum of the intensity of all 
reactions, ao(x), gets large: since all the jump times have to be included in the time 
discretization the corresponding computational work may become unaffordable. 
Indeed, we have that the mean value of the number of jumps on the interval (t, t+r) 
is approximately ao(A t )r + o (r). 

Assumption 2.1 (Bounded population). We here assume that species can only 
be transformed into other species or be consumed. This means that the state X t G 
Z + will be bounded from above by a hyperplane intersecting the point X$ and the 
coordinate axes, i.e. 

X t £ n(A , n) := {x^7L%:n- (x - X ) < 0}. 

for some vector n £ with strictly positive coordinates, i.e. n > 0. 
In other words, all the stoichiometric vectors should satisfy 

n-v,< 0, j = 1, . . . ,M. 
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Remark 2.2 (Preventing negative populations). To prevent the state from becom 

2 j 



ing negative after a jump we must have that aj(x) — for x + Vj ^ lA. For 



the same reason we also have that aj(0) = 0. For common chemical reactions the 
propensity functions can be modeled as polynomials of the form 

d | 

which satisfies the above criteria, see [2, 11]. Here, := min(2/y,0). From the 
expression (2.3) it also follows that the propensity is monotone in l\, i.e. the 
gradient of the polynomial function dj is non-negative. 

Remark 2.3 (Boundary values). From Remark 2.2 it follows that Equation (2.1) 
has a natural boundary condition at the boundaries where any component of the x 
vector is zero. 

3. The tau-leap method 

To avoid the computational drawback of the SSA, i.e. when many reactions occur 
during a short time interval, the tau-leap method was proposed in [10]: Given a 
population X t , and a time step r > 0, the population at time t + r is generated by 

M 

(3.1) X t+T = X t + J2 Vj-PM(X t )T), 

3=1 

where V{aj{Xt)r) is a sample value from the Poisson distribution with parameter 
aj(X t )r, indicating the sampled number of times that the reaction j fires during 
the (t, t + t) interval. 

The tau-leap method (3.1) is based on the observation that if the propensity 
is constant between t and t + r, the firing probability in one reaction channel is 
independent of the other reaction channels. The total number of firings in each 
channel is then a Poisson distributed stochastic variable depending only on the 
initial population X t . Also, Equation (3.1) is nothing else than a forward Eulcr 
discretization of the SDE corresponding to the Kolmogorov backward equation 

(2.1) , (or the chemical master equation), driven by the Poisson random measure, 
cf.[}- 

In the following we let X t denote the exact process and X t the tau-leap approx- 
imation. 

Remark 3.1. Although the path generated by (3.1) seems to only be defined for 
discrete time steps, it can be extended to the continuous time interval. Indeed, any 
intermediate time steps in the (t, t + r) interval can be defined as a Poisson bridge 
with fixed endpoints X t , X t + T and fixed intensity a~j := aj(X t ). 

After freezing the propensity to a{X t ) in each time step [t, t + r], the error in the 
tau-leap method (3.1) comes from the variation of a(X s ) for s e (t,t + r), where 
X s is the true process starting at X t — X t . To take this into account it was in [10] 
proposed choosing local time steps following the leap condition 

(3.2) |A aj (A t )| := \a.j{X t+T ) - a 3 {X t )\ < ea (X t ), 

for a given control parameter < e <C 1. In order to avoid unnecessary sampling of 
the Poisson random variable a pre-leap check can be done by first doing a Taylor 
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expansion 

A aj (X t ) = a 3 [X t +Y,ViPi(oi(Xt)T)) - a 3 (X t ) 



i=l 
M 



« Va i (X t )-^i/ i P i (a i (X t )r), 

i=l 

and approximating the mean and variance of Aa.,(a;) by 
1 M 

-E[Aa,-(a;)] « y^qj(x) (Vflj-(aj) • vC) =: fij(x), j = l,.. .,M, 



and 



— Var[Aaj(x)] ~ y^a,-(a;) (Va,j(x) ■ Vi) 2 =: ct^ie) 2 , j = 1, 
T i=l 

respectively, see [13]. The leap size is then chosen as 

'eao(Xi) e 2 a (X t f 



,Af, 



(3.3) 



mm mm 



i=l,..,M V 1^(^)1 »iW 



which implies that |E[Ao,-]| < ea Q and VarfAeij] < e 2 al, so that the leap condition 
(3.2) holds in some statistical sense. The leap size control (3.3) has the advantage 
that it can be computed before each step is taken; however, it relies on the ad hoc 
parameter e and only controls the local error. 

Our goal is here to develop a rigorous a posteriori global weak error estimate for 
(3.1), with computable leading order term, that can be used to control the error in 
a more systematic way than the leap condition (3.2), and that can in future work 
be used for efficient adaptive time stepping algorithms. In Section 3.1 we present 
a procedure that eliminates the possibility of negative populations X t , similarly to 
[1], and in Section 4 we develop error estimates that are finally tested numerically 
in Section 5. 



3.1. Avoiding negative population. The regular tau-leap approximation (3.1) 
suffers from the undesirable property that X t can become negative. To prevent such 
unphysical behavior, which is a result of the approximation and not the process 
itself, the idea is to adaptively adjust the time step r to avoid negative populations 
and at the same time leave the distribution of X t unchanged. Similarly to [1], this 
is done by post-leap checks such that if any component of X t + T becomes negative, 
a new step X t+T / 2 will be sampled using a conditional Poisson distribution, i.e. 
a Poisson bridge. To describe this procedure, we first introduce independent unit 
rate Poisson processes Yj(-), and define their corresponding internal times as 

A/ : / a,:\.u\,. j 1 M. 

Jo 

The state of the chemical system at time t then satisfies 

M 

X t = X + J2"jYj(4), 
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see [1], i.e. each reaction channel is described by a unit rate Poisson process with 
an internal time given by X{ . From the definition of Yj we have Poisson distributed 
increments 

Y j(H+t) ~ YAK) = Yi U + J^ +T a 3 (X s ) ds) - Y 3 (Xl) 

= T>i (J T aj(Xs)ds\ , 
and since we have for each tau-leap step that 

[t+T 

(3.4) J aj (X,) ds = a 3 {X t ) T := AX{, 

the tau-leap method can be written in terms of increments in Yj, i.e. 

M 

(3.5) X t+T =X t + J2 <^V;,:AA/i. 

3=1 

where 

(3.6) AYj(-) :=V 3 {-). 

As the tau-leap method steps forward in time we build up a history of samples 
for the driving process Y, {(A^, Y^)}™1 , such that = Aq < X\ < . . . < X J n ., and 
= Yq <Yf < ... < Yji., by summing increments (AX J t , AYj (AX J t )) . Starting at 
X t and at the last value (Ai,Y^) in the history, a regular tau-leap step samples 
an increment AY t 3 := AYj(Axf) and saves (X j k + AX{,Yg + AY t j ) to the history. 
Let B(n,p) denote the binomial distribution. As soon as a negative population 
Xt+T < is encountered, a new step X t+T / 2 is calculated by the Poisson bridge 

M 

X t+ r/2 = X t + J2 VjB(&Y t j , 0.5), 

3=1 

and (XI + 0.5AX{,Yj! + B(AY t j ,0.5)) is added to the saved history. If X t+r/2 is 
non-negative we move forward along the £-axis, otherwise the step is halved once 
again. When traversing the physical time axis t we have to check if there already 
exist samples to the right on the internal time axis A, e.g. given the position 
A^. and a future value A^ , 1 the physical time step r must be adjusted such that 
we end up on X k , x for at least one j and to the left of X J k+1 for the remaining 
j. For the reaction j that ends up on X J k+1 we use the corresponding Y k+1 , and 

for the other reactions, a bridge between (X k ,Y k ) and (A^ +1 , Y^ +1 ) is sampled. 
Algorithm 3.1 describes in detail one step of the Poisson bridge tau-leap method 
in the interval (t, t + r). To speed up the algorithm we approximate, for large np, 
the binomial distribution B(n,p) with the normal distribution J\f(np,np(l — p)) 
rounded to integers and multiplied by the indicator function with support in [0, n]. 
We here apply this approximation whenever np > 10 4 . 

Remark 3.2. The post-leap check performed by the Poisson bridge tau-leap method 
in Algorithm 3.1 will guarantee non-negative sampled steps. There is however still a 
probability that components of the continuous tau-leap process may become negative 



8 



J. KARLSSON AND It. TEMPONE 



Algorithm 3.1: One step for the Poisson Bridge Tau-Leap method. Unless 
stated otherwise, it is assumed that steps in the algorithm containing the reac- 
tion index j are performed for all reactions j = 1, . . . , M. 

Input : Initial values (t,X) and (X*,Y*). Coefficients Vj and dj(x), step 

size t, and an ordered history of samples Aj := {(Xj, Y/)}j=o) where 
<> K • Ai • ... A; and = Y* < Y( < . . . < Y*. . 

Output: An ordered history of samples Aj, the current value of values for 
(X*,Y*), and the steps {(tj, Xj)}?l taken in the interval [t,t + r]. 

Set local time i/ oc = 0. 
while (tioc < t) do 

Let aj := a.j(X), and let kj be the index of current X J in the sample 
history Aj. Set local time step t; oc := r — £/ oc and corresponding internal 
time step AA J := a~jTi oc . 
for (j = 1, ...,M) do 

Save a temporary local time step t? := t; oc for reaction j. 
if (kj is last index in Aj) then 

Sample a new Poisson value AY' J := ^(AA- 7 ), and add 
(A^. + AX* ,Y£. + AY 3 ') to the history Aj. 
else 

if :A; ■ A A' AJ . : i then 

Set local time step for reaction ?': rf := 3+ _ — - x . 

" ^ toe aj 

Set AA< : - X{. and AY' := - Y*.. 

else 

Sample a Poisson bridge value 

AY 3 := B(Y/ +1 - Y"/., j AA \ i ) from the binomial 

J J fc j + 1 k j 

distribution, and add (X{. + AX j , Y£ + AY 3 ) to the history Aj. 

end 
end 
end 

Set local time step to be the minimum physical time to next sampled 
value in history: r/ oc := mim, Tj oc . 

while (any component of X + J^j VjAY ] is negative) do 

Divide step size by two, i.e. ti oc := 0.5t; oc , and AA J := 0.5AA- 7 . 
Sample a bridge value Ay' := B(AY* , 0.5). 
Add (x{. + AA-'. V,; + AY*) to Aj. 

end 

Accept new value X <— X + J2j VjAY* , increase internal time 
A J <— Aj.. + AA J , and increase local time ti oc <— ti oc + n oc . 
end 



at points between the sampled steps, as discussed in Remark 3.1. To limit this effect 
we can introduce a pre-leap check that adjusts the time step r such that 

(3.7) P(A > t ( ; ) T < | X t ) < e, i=l,...,d, 
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for some small e > 0, see [1.5] where the exit probability (3.7) is approximated by a 
normal approximation of the tau-leap step (3.1), leading to a quadratic inequality 
in t. 



4. Accuracy and error estimation 

To develop computable global error estimates for the Poisson bridge tau-leap 
method we start with the following (non-computable) error representation based 
on the backward Kolmogorov equation (2.1), along the continuous-time approxi- 
mate tau-leap paths X t , defined in (3.1) and Remark 3.1, and the difference in 
propensities: 

Lemma 4.1 (Error representation for the tau-leap method). Assume that aj and 
u are defined for x £ Z d and that u solves the backward Kolmogorov equation (2.1) 
in Zi x [0, T] . For deterministic time steps we then have 



M 



E[g(X T )-g(X T )] = £ E 

3=1 



(a-j - Oj)(Xt) (u(X t + vj,t) - u(X t , i)) dt 



E 



4>{X u t)l {Xt& d _} dt 



where 



M 

<t>(X t ,t) := d t u(X t ,t) + aj{X t ) (u(X t + Uj,t) ~ u(X t , tfj , 

3=1 

and Z_ := % d \ lA is the set where at least one component is negative. 

Proof. Since Xq = Xq we have 

E [g(X T ) - g(X T )] = E [u(X , 0) - u(X T , T)] 

(4.1) 



= E 



-(d t + Cx)u(X t ,t) dt 



where the "infinitesimal generator" for X t is defined as 

M 

C x u{X t ,t) \=^aj{Xt){u{X t + Vj,t)-u{X u t)). 



3=1 



Note the abuse of notion here: dj, and correspondingly C x u, is from Remark 3.1 
only defined along paths. Adding and subtracting 



A I 



C x u(X t ,t) :=Y,a {X t ){u{X t + Vj ,t) -u(X t ,t)), 

3 = 1 



gives 



E [g(X T ) - g(X T )] = E 



{Lx - C x )u{X u t) - 4>{X u t)dt 



From the Kolmogorov backward equation (2.1) we have <fi — for X t £ Z^f. which 
gives Lemma 4.1. 

□ 
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4.1. Relative error. Assume that we have x £ Zl particles and propensities 
a,- (a;). Introduce the scaling z :— x/-f where the scaling factor 7 > is related to 
the initial number of particles Xq, e.g. by taking the Euclidean norm 7 := ||Ao||. 
The relative error can now be defined as 

(4.2) e:=E[g(Z T )-g(Z T )], 

where the scaled Z t -process is represented by the Kolmogorov backward equation 

M 

d t u + ^2aj(z) (u(z + Pj,t) -u(z,t)j = 0, in Z f {/ 7 x [0,T), 

(4.3) i=i " ^ ' 

u = g, on Z^/ 7 x {T}, 
and approximated by the rescaled tau-leap method 

M 

(4.4) Z t+T =Z t + J2 *iPj (2j (Zt)r), 

i=i 

with the scaled propensities cij(z) := aj(x) — 0^(72) and jumps Vj :— 1^/7. 

Example 4.2. Let g(x) = x r := H* =l and \r\ := ^ d =1 \n\ > 1 w/iere r g Z^_ 
is a d- dimensional multi-index. Then 

E[g(X T ) - g(X T )} 



9(Xo) 



E[g(Z T ) - g(Z T )] = 
for some min^ X^ < 7 < maxj Xq . 

Assumption 4.3. Following Assumption 2.1 and Remark 2.2 we here assume that: 

(1) The process Z t £ Z^/7 is bounded from above by a hyperplane passing 
through the point Zq and with normal n > 0, i.e. 

Z t £ il{Z , n) := n( 7 A , n)ft = {z £ Z^/ 7 :n-(z-Z )< 0}. 

Consequently, the tau-leap process takes values in the larger set 

Z t £ fl(Z , n) := {z £ Z d /7 :n-(z-Z )< 0}. 

(2) The propensity aj : Z^ — > M is non-negative and non- decreasing in each 
component. Also, we have that dj(0) — 0, and aj(x) = if x + v } ■ (£. Zi. 

The goal is now to show an a priori bound for the discretization error e that 
is independent of the scaling factor 7, see Theorem 4.4. The proof is based on a 
Taylor expansion and bounds on the value function u and its derivatives. 

Theorem 4.4 (Bound on relative error). Let the process Z t and the propensities aj 
satisfy Assumption 4-3. Assume polynomial propensity functions of order \pj\, with 
Pj £ Zl, and let g £ C 2 (R d ). Also, assume that the value function u given by (4.3) 
and its first two spacial derivatives are bounded in Tl(Zo,n) x [0, T], independently 
of the scaling factor 7. 

For a time step t — hj~ s , with h > and S = m&Xj(2\pj \ — 2), the relative error 
(4.2) for the Poisson bridge tau-leap method is then bounded by 

\e\ < Ch, 

for some sufficiently large 7 and where the constant C > is independent of h and 
the scaling factor 7. 
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Proof. Throughout the proof C will denote a non-negative constant value, not 
necessarily the same at each instance. To show the theorem we note that from the 
error representation in Lemma 4.1, the relative error can be expressed as 

M ftn+l 



(4.5) 



where 



-EE 

j = l n=0 



E 



(fflj(Zt) - a,j(Z tn ))Dju(Z t ,t) | Z tn 



dt 



w-i ,t„ +1 



E 

n=0 



E 



tn 



(A n u(Z t ,t)) \ { z t & d _/i} I Z t„ 



dt. 



A I 



A„ :=d t + Y,~ a i&n) D r 
For the first term in (4.5), Taylor expanding the propensity a around Z tn gives 



\Pi\ ?,( a ^(7 \ 

/ .- — — E 



(Z t - Z tn ) a iD 3 u{Z u t) 



where <x,- € Z+. 

According to Lemma 4.5 the value function and its two first spacial derivatives 
can be continuously extended to M. d x [0, T]. The mean value theorem then yields 



{Z t - Z tn ) a i (u(Z t + v jt t) - u{Z u t)) 



1 



(Z t - Z t X 3 [ d z u{Zut)v, + d 2 z u(^t)i) 3 







) 





[Z t - Z tn ) a id z u{Z tn ,s)i> 3 z.,,. 
E \{Z t - Z tn ) a3 (Zt - Z tn )d*u(&, t)vj 
E uT(Z t - Z t „) a >tfu(tut)u4 Z t , 



vJ{Z t ~Z t X 3 d 2 z u{^t)v 
ii i fin 



p + + 7 



for & = aZt + (1 - a) (Z t + £>,,•) for some a £ [0, 1] and £ 2 = fiZt n + i 1 ~ P) Z t for 
some f3 G [0, 1]. 

From the tau-leap approximation of the process Z t we have 

A I 



E 



= Y.v[vdHZu){t-t n )t A W 

3=1 

AI I ay I 
M I ay I 



3 
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for some positive constants Cp j that do not depend on 7. Choosing t = /17 A for 
/i, (5 > 0, and without loss of generality assuming that the propensity is a monomial 
aj (x) = Cx p i with C > 0, pj £ Z| and |^| > 1, then 



a i (Z t jftrft|i/f| = a J -(7^J 



= a i (7^ t jft^|^ a3 '|7-^-l^l 



Using the continuous extension in Remark 4.6 and assuming the value function u 
and its first two derivatives are bounded independently of 7 inside H(Zq, n) x [0, T], 
i.e. neglecting the probability terms in the growth condition in Lemma 4.11, leads 
to continuity and boundedness in the extended domain {z £ R d : n ■ (z — Z ) < 0}. 
Together with 6 — maxj(2\pj\ — 2) and h < 1, this gives 



\Pj\ rX a i)(7 \ 

a A 



lay 1=1 



\Pj\ K'l 

< c V V 7 lftHH^+&(lftl-<5)-i/j& 

K|=i&=i 

< 7 max i(|Pil)-2+/3j(max J (|p : ,-|)-(5)^ 
1^1=1 ft-=l 

< Q 7 max j(|Pjl)- 2 +(2-max j (|p ;j |))^ 

1^1=1^=1 

< Ch, 



for Z tn e IL(Z ,n). 

Performing the same analysis on the error contribution from the terms 1^ n and 
T^ 11 n gives the estimate in Theorem 4.4. Note that the additional (Z t — Z tn ) term 
in^ J)fl and the extra v s in PJj n yields that P a ^ n = O (1), R] jin = O (7- 1 ) and 
T^J n — O (7 -1 ), so for large 7 the first term in the Taylor expansion is dominant. 

For the second term in (4.5) we have that since the value function and its deriva- 
tives are continuously extended and bounded in II(Zo,n) 



A I 



A I 



A n u(Z t ,t) = [ dt + ^ajiZtjDj u{Z u t) < 7 _1 X^a^ZtJ- 
3=1 / 3=1 



and thus 



A I 



E 



A n u(Z t , i)l{ ZtGZ d / 7 j I Z tn 



< 7" E 



L {z t ezl/ 7 } 



j=i 
,1/ 



(4.6) 



X P (Z t e Zi/ 7 I 4)E C i a i^«) 



2=1 



71/ 



< 7- 1 P (z« < I Z tB ) C,a,(Z t J. 

i=l 2=1 



COMPUTABLE WEAK ERROR EXPANSION FOR THE TAU-LEAP METHOD 



13 



Neglecting positive jumps we obtain 

M 

Z t W = Z«+X;^i(5i(^)(*-*n)) 
3=1 

> Zg + v^V{a {Z tn )T) =: z\ l \ t e {t n ,t n+1 ), 

where i>W := min :) {min(i>j i \ 0)} < 0, and a : = Ejli For > the proba- 
bility of negative populations in (4.6) can then be approximated by 

P (z^ < I Z t \ < P (z ( / } < I Z tn 



< P 



E 



(4.7) fe=W 



/ 

P(5o(Z tn )r) > -Z®/u® 
, v . ' 

V =:g>0 

fc! 



E 



T 



(5 (Z t Jr) fc+ ^ 
ffc -I- \a\ 

k=0 



-a {Z tn )^ 



< ^ ^ \q\\' ' /c! ' 6 
fc=0 1 ^ 1 

(5 (Z t Jr)^ 

(C 7 )l ' 

where [•] denotes the ceiling and C is independent of 7. For Z^ = 0, the proba- 
bility P {z[ i] < I Z t \ is zero since a 3 -(Z t J = if ^ < 0. 

Combining (4.6) and (4.7), and using o (x) < Cx ma ^ M and r = /if- 2 ™' 1^1 
gives 

E A„u(Z t ,t)l { ^ eZ d /7} <Cj~ 1 a o {Z t J^-0^ — 
(4-8) < C y-^l ft l-i(^r^ 



,c 7 



(C 7 )! 



so the remainder term (4.6) can be neglected. 



□ 



In the theorem above we used that the relative value function u defined by (4.3) 
and its derivatives can be continuously extended to R d x [0, T], and are bounded 
independently of 7 in the domain Yl(Zo,n) x [0, T]. In the remaining part of this 
section we will motivate why those assumptions make sense. 

Lemma 4.5 (Extension of the value function and its derivatives). Assume that the 
process Z t and the propensities aj satisfy Assumption 4-3, and that aj and g are C 2 
on M. d . The solution u to the Kolmogorov backward equation (4.3) is then locally 



14 J. KARLSSON AND R. TEMPONE 

bounded and smooth in time on Z^/7 x [0,T] and can be continuously extended to 
R d x [0, T]. Also, the derivatives d z ii and d 2 u are locally bounded and smooth in 
time on the lattice Z^/7 and can be continuously extended to M. d x [0, T]. 

Proof. First, we continuously extend dj such that it is non-negative, monotone in 
each component and C 2 on K d , see e.g. Example 4.7 and 4.8 below. Given such an 
extension of the propensity, the value function u can be continuously extended to 
the domain M. d x [0, T] by solving the Kolmogorov backward equation (4.3) on the 
shifted lattice (Z d + + e)/ 7 x [0,T], e £ R d . 

Note that, for any point Z t £ Z^/7 we know that Z s € Il(Z t ,n) for s > t, 
and (4.3) can be understood as a linear constant coefficient system of ordinary 
differential equations on a finite lattice, with a unique bounded smooth solution it. 

The derivatives v := d z u £ M. d and w :— d z u £ M. dxd can be defined on Z+/7 x 
[0, T] and satisfy the equations 

M M 

d t v + a j (z)D j v = -^2 ^a j (z)D j u, in Z d + /j x [0, T), 

3=1 3=1 

v = V.g, in Z^/ 7 x {T}, 



M 

dtw + dj(z)Dj€> = 
3=1 

M 

^2(-2Vaj(z) ■ D 3 v - V 2 a 3 (z)Dju), in Z^/7 x [0,T), 
3=1 

w = V 2 g, in Z%h x {T}. 

Using that g £ C 2 (M. d ) and that u is locally bounded and smooth, we can by the 
same argument as for u above see that v and w are locally bounded and smooth in 
time. 

By shifting the grid in the same manner as was done for u, the solutions v and 
w can be continuously extended to K d x [0, T]. 

□ 

Remark 4.6. In addition to the above lemma we can also continuously extend the 
value function u(z, ■) to real negative values z £ Ml such that u, d z u and d 2 u vanish 
for Zi < ~l < 0. This extension is used for the Taylor expansion in Lemma 4-4- 

Without loss of generality consider the extension to the domain M^, with sub- 
domains A := {z £ (—1,0) x (0,oo)} 7 B := {z £ (-L0) 2 } and C :— {z £ 
(0,oo)x(-f,0)}. 

First let u(z, ■) = on z £ M_ \ {A UBU C}. Next, extend the value function 
to the domain A by third degree polynomials in the z\ direction, for all Z2 > 0, 
matching u, d Zl u and d Zl u on z\ = —L and z\ = 0. Similarly, fit polynomials in 
the Z2 direction in domain C. In the domain B a C 2 -extension can be made by 
transfinite interpolation using it, d z and d 2 u on dB, see e.g. [3, 14, 27]. 



and 



(4.10) 
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Example 4.7. A natural extension to the propensity function dj(x) := x, x G Z, 
is 

{0, x < 0, 

Sj (ac), < x < 1, 
3C, X > 1, 

/or a; € R, where dj is a positive and monotone C 2 function chosen such that 
dj e C 2 (R), e.g. aj(x) := 6a; 3 - 8x 4 + 3x 5 , see Figure 4.1. 

Example 4.8. The propensity function a(x) :— x(x — 1) is negative in (0, 1). A 
natural extension is thus to let 

(0, x<l, 
a(x) := < a(x), 1 < x < 2, 

[x(x-l), x>2, 

for x G M, where &j is a positive and monotone C 2 function chosen such that 
a,j G C 2 (R), e.g. aj-(x) := 9(ac - l) 3 - 11 (as - l) 4 + i(x - 1) 5 ; see Figure 4.1. 




(a) Example 4.7 



(b) Example 4.8 



Figure 4.1. Extended propensity functions in Example 4.7 and 4.8. 

In Lemma 4.5 we saw that u and its first two derivatives are locally bounded 
in R d x [0, T], however, this bound depends on 7. The next step is to show under 
what conditions we can find a bound that does not depend on 7. First, we start by 
generalizing the polynomial form of the propensity (2.3) in Remark 2.2 into bounds: 

Assumption 4.9. Given some multi-index pj G and positive constants < 
C l L < C\j, < C' L < C'jj and < C'l < Cy, assume that the propensity has the 
bounds 



(4.11) 



for i = l,. 



C L (x + Q)7 < <Cu(x + 0)? , 

c L (x + Q) p 4~ ei < d Xiaj (x) <c u (x + Q)T ei > 

Cl(x + < d Xzk a 3 (x) < C'iix + QT^^ 

.,d, j — l,...,M and x G Z 1 ^. Here, (-) + := max(-,0), Qi 



(v,ij + l)l^. j< o, and ej indicate the unit basis vectors in 
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In the next step, we introduce a weight definition that in Lemma 4.11 will allow 
us to do find weighted estimates for the value function and its derivatives. 

Remark 4.10 (Weight definition). From the Assumption 4-9 on boundedness of 
Z t , there exists a vector n > in K d such that i/j ■ n < for j = 1, . . . , M. Take a 
weight function y(z) :— ip(n-z) > where ip '■ (t m - m , +oo) — > K + is a smooth, strictly 
decreasing function. Observe that we have by construction y(z + Vj) — y(z) > 0, 
for all z £ II(Zo,n). In particular, if g has polynomial growth g{z) < C(z ■ n) r °(z ■ 
n + l) ri we can have the product yg uniformly bounded, for instance by taking 
4>(t) = (t + t )~ r "(t + l)~ ri , with some -£ min = <o > and r ,ri > 0. 

Now we are ready to present weighted estimates that are uniform with respect 
to the scaling parameter 7, and on which Theorem 4.4 for the a priori relative error 
bound is based. 

Lemma 4.11 (Growth condition). Assume that the assumptions on boundedness 
of Z t and the conditions on aj in Assumption 4.3 hold. Also, assume that the 
bounds (4.11) on the propensity and its derivatives in Assumption 4- 9 ore satisfied. 
Let \p\ be the maximum reaction order in the system, i.e. \p\ := max^ \pj\ for the 
multi-indices in (4.11), and assume that for each of the components Xi, there exist 
a single component reaction whose propensity has order \p\. Finally, assume that g 
is C 2 on R d , non-negative with g(0) — 0, and that for the weight functions 

y(z) := z ■ n H — max ■ nil (l + z-n) 71 , 
V 7 ) 

yi(z) := y{z) ( z ■ n H — max \uj 
\ 7 J 

y2(z) :— yi(z) ( z ■ n H — max \va ■ n 
V 7 i 

the quantities \yg\, |yiVg| and |j/2V 2 g| are bounded. 
We then have the bounds 

\u(z,t)\< ( 



2/0) 

(4.12) \d z u{z,t)\< -^- + -^ r -(^z-n)^p(- <Z T -n<r] Z t = z] 
yi{z) yi(z) V7 / 



\dlu(z,t)\ < -^ + ^L( lz . n f\v\ V ( e _ <ZT . n<71 
y2{z) 2/2(2) V7 



Z t = z 



for (z,t) € LI(Zo,n) x [0,T] , for some 6,rj > and positive constants C, C and 
C" that are independent of"/. 

Proof. Throughout the proof C will denote a constant value, not necessarily the 
same at each instance. Also, for simplicity we will assume that the propensity 
has the form aj(x) :— Cx Pj . In the general case with bounds (4.11), the proof is 
essentially the same. The proof is divided into three sections corresponding to the 
bounds in (4.12). 

Bound on u: Let U(z,t) := y(z)u(z,t). Since g has polynomial growth we can 
take y as in Remark 4.10 with 

to = max \vj ■ n\ :— \D* ■ n\. 

j 
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Also, since g(0) = we have rp > 0. Then we have 



M 

d t U + o-j (yu(z + pj,t) - Uj = 0, in fl(Z , n) x [0, T), 



3=1 



where \\U(-,T)\\ 



subtracting Y^jLi ®j( z )U{z + Vj, t) yields 



\yu{-,T)\\, 



U = yg, on fl(Z , n) x {T}, 
\V9\\oo by definition is bounded. Adding and 



M 



M 



(4.13) 



dtU + 5,jDjU = djii(z + Oj,t)Djy, 

3=1 3 = 1 

with the stochastic representation 

T m 

U(z,t) = E y(Z T )g(Z T )~l ^ a j {Z s )u{Z a + %, s)D j y(Z a 

Jt j= i 

Since the right hand side of (4.13) is non-negative we obtain 

l|C(-,*)Iloo < ||cr(-,r)||oo < c, 

and therefore 



Z t = zds 



u(z,t) < Cy(z)- L < C [z-n + 



7 



(l + z-n) r 



which is the first estimate in (4.12). 

Bound on du: For the bound on the first derivative v := d z u € R d , we consider 
the weighted function V(z,t) := yi(z)v(z,t) £ M. d with the weight 



yi( z ) : = y( z ) [z-n + 



\v*-n\ 



7 



z ■ n + 



\v* ■ n\ 



7 



l-ro 



(1 + z ■ n)~ 



such that the product yiV g is bounded uniformly in z and 7. Define an auxiliary 
function 

M 

R(z, t) = a,j(z)v(z + Uj,t)Djyi(z). 
3=1 

This yields 

M M 

d t V + ^2 <ijDjV = R(z, t) — ^2 Vaj(z)yi(z)Dju(z,t), mfl{Z ,n)x [0,T), 
3=1 3=1 

V = yi Vg, on Ii(Z , n) x {T}. 

Denote 

Pj( z ) ■=a j (z)D j y 1 (z)y 1 (z + £>j) - \ 
(4.14) *£, 

3=1 

and observe that /3j(0) = and Pj >Q. Since 

M 

Ji(z,t) = 5^V(« + i? i ,t)/8i(«), 
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we have 

M M 

d t V + ^(dj - Pj)DjV - f3V = -yi ^VajDjU, inIT(Z ,n) x [0,T), 

V = yi Vg, onn(4n)x{T}. 

The solution V can now be stochastically represented as 
(4.15) 



V(z,t) = E yi (Z T )Vg(Z T )e-f* ^ ds 



Z t = z 



=Vi 



E 



M 



yi (Z s )e-f'^ dt ' ^ajiZ^DjuiZ^s)^ 

3 = 1 



Z t = z 



=v 2 

where Z t is a modified jump process with smaller propensities 

aj{z) ■= (a,j-Pj)(z) = aj(z)(l - D j y x {z)y l {z + Oj)' 1 ) = a j {z)y 1 (z)y 1 (z + v^ 1 . 

Indeed, since y\(z + > yx(z) > by definition, we have < cij(z) < a,j(z). 

Observe that since the product \yiVg\ is bounded, we have \Vi\ < C. We now 
consider V2, and first focus on a lower bound for the functions /3j defined in (4.14). 
From the assumption aj{x) := Cx Pj we obtain 



f3 J (z) = Cj^zP^ J (z), 



and 



£j(z-n) := D jyi (z)yi(z + i>j) 
= 1 -11- |f>j " n| 



z ■ n + \v * ■ n\ 



ro — l 



1 - 



1 + z ■ n 



Observe that £j is strictly decreasing with respect to the product z ■ n on (0, +00), 
and that £(+00) = + . Since the values z G H(Z Q ,n) and v are bounded by 
assumption 2.1, we conclude that £j(z) > £0 > and then 

for all z £ n(Z ,n). 

Given r\ > 0, let us now introduce a family of time intervals, 



Let 



/* :={se[i,T] :Z s -n>r]}, 



\p\ = max \pj I 



and observe that, thanks to the assumption on the existence of all the single com- 
ponent reactions with order \p\, 



(4.16) 



/3(Z fl )da>6,C|J'|(7»7) 



\v\ 
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On the other hand, we have 



\®zi aj {z)y\ (z)Djii(z, s) | <C\d Zi aj (72)2/1 (z)y(z) 



(4.17) 



= C 7 IP 3 'I; 



•Pj- 



z ■ n + 



<C^ p t\z-ri)\ p i\- 1 [z-n + 



\v* -n\ 

1 



Let < k < 1 denote a time fraction. Then thanks to (4.16) and (4.17) the term 
V 2 in (4-15) can be estimated as 



\V 2 (z,t)\ < C 



E 



7 l p l(i 5 -n)l p l e -«oC(7r,) lpl |^l 



Z t = z 



ds 



< (7 7 bl / e -ZoC ( 7 »?) |p| «(r-*)ds 



(Z s • n) |p| l|/s| <K(T _ s) ds 



C 7 |p| E 



Z t = z 



1 _ e -CoC ( 7 ,,)l p l K (T-t) /.T 

<C goC ^ |p|K + C7 |p| J P(l^l < «(T - s)\Z t = z)ds 

Observe, that since the process Z s ■ n is non increasing, we have 

i^K^T-S!) < 1 |/^|<«(T- S2 ) < 1 z T -n< V ' for s i < s 2 < T, 
and then, for a given constant 8 > 1, 



7 



bl E 



(Z s • n) |p| l|/s| <re(T _ s) 



< 7 i p iE 



< 6» |p| E 



(Z s -n)\Ph 



(4.18) 



1 Z T -n<« ds 



+ 7 N( z . n )blE 



ds Z t = z 

Z T -n<n \}-Z B n< 
Z t — z 

ds 



L Z s -n>i 



ds 



Z t = z 



- <Z T -n<7] 



Z t = z 



< {T-t)6 M P (z T -n< 
+ (T-t) 1 M (z-n) M P 



Z, = z 



< Zt ■ n <r) 



Z t = z 



Observe that when the process Z s ■ n is within a distance 0(1/7) from zero 
there are only bounded contributions to the derivative of u. Therefore, what really 
deteriorates the derivative estimate is the time spent between 0/ 7 and 77. 

Bound on d 2 u: Similar to the bound on du we consider the weighted function 
W(z,t) := y 2 {z)w{z,t) = y 2 {z)d 2 z u{z,t) € R dxd with the weight 



y 2 (z) = Vi{z) [z ■ n + 



\v* ■ n\ 
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such that the product y 2 V 2 <? is bounded uniformly in z and 7. Following the same 
procedure as for V we obtain from (4.10) that 



d t W + ^{aj - Pj)DjW - PjW = -R, in IL(Z , n) x [0, T), 

VF = y 2 V 2 .g, on n(Z , n) x {T}, 



where 



A/ 



(4.19) /3(2) /3,(z):=a J (z) J D, 2/2 (z)y 2 (z + ^)- 1 , 

and 

M 

i?(z, t) = ^ 2y 2 (z)Va 3 (z) (8) t) + V 2 a.,(;z)y 2 (z).D,-u(2:, *), 

j=i 

with '(g)' denoting the tensor product. This gives the stochastic representation 
formula 



W(z,t) = E y 2 (Z T )V 2 g(Z T )e-^ 



/ t T /3(z s )d s 



Z+ = z 



(4.20) 



+ E 



y2(^ s ) e" 



=:Wi 



-/ s T /3(Z t ,)dt' 



R(Z S , s)ds 



Z t = z 



:W 2 



where the modified process Z t has propensity a,j(z) := cij(z) — j3j(z). 

Since \y 2 V 2 g\ is bounded we have that \Wi\ < C. For W 2 we note that (4.16) 
holds also for the new f3 in (4.19) and the new modified process Z t . Also, we have 
the bounds 

\d Zi aj(z)y 2 (z)Djv(z,s)\ < C\d Zi aj{jz)y 2 (z)\ 

■ 7 |p| (z • n) |p| P (o^ 1 < Z T -n<ri\Z s = z) yi(z)" 1 



C 7 M z Pj-e< ( z -n + 



\v* ■ n\ 

7 

7 bl( z . n )blp ^ 7 -i < z T -n<r)\Z a = z^ yi(z) -1 
)2|pj|p < Z T -n <r]\ Z s = z) , 



and similarly 

l^ Zfe aj( 2; )y2(z) J D J u(z,s)| < C\d ZtZk a ] {^z)y 2 {z)/y{z)\ 



= C7^"l; 



,pj-ei— ek 



z ■ n + 



1 1/* •n| 



< C 7 ^l(z-n)lwl- 2 ^.n+^-^y 

< C7l p 'l(«-n)'wl. 
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In the same manner as in (4.18) we obtain 



\W 2 (z,t)\ < C 



7 |p| (Z s -n)l p le 



Z t = z 



C J E 7 2|wl (Z s • n) 2| ^!p (Vr 1 < Z T ■ n < r) \ Z s , = Z t 



x e 



-?oC( 7t? )' p l |J*| 



Z t = z 



As 



< C + C7 |p| (z • n) M P (e-/- 1 <Z T ■n<r)\Z s = 



+ C~f 2lpl (z-n)~ zm P [9-f- 1 < Z T ■ n < i] \ Z s = z 



and in particular if the quantity (4.21) is uniformly bounded in 7 we have that 
\W\<C. 

□ 



Corollary 4.12. Assume, in addition to the assumptions made in Lemma 4-11, 
that given z > there exists rj(z) s.t. for a sufficiently small constant 9 > the 
quantity 



(4.21) 

is uniformly bounded in 7. Then 



- i' (: • ») '' P ( ^ <Z T -n<r) 



Z t = z) <C 



\ u ( z , t)\ < C [ z ■ n + max \vj ■ n\ ] (1 + z • nf 



\d z ii(z, t)\ < C ( z ■ n + max |z>j • n| 



r — 1 



(l + z-n) r 



d z u(z, t)\ < C [ z ■ n + max |£j • n 



1-0-2 



(1 + 2 • n) T 



/or (2,t)en(Z„,n)x [0,T]. 

Tfte assumption (4.21) is trivially true whenever n-Vj =0 for all j = 1, . . . , M, 
e.g. /or a system with only reversible reactions of the form X\ ^ X2. In that case 
the product Z t ■ n — z ■ n remains constant for all t < s < T. 

4.2. Computational work. The results in the previous section now give us the 
possibility to judge in which regimes and for which propensities the Posson bridge 
tau-leap method is expected to be more efficient than the SSA. As mentioned in 
Section 1, the computational work of the SSA is roughly inversely proportional to 
the total propensity, and becomes intractable as the number of particles grow. The 
tau-leap method, on the other hand, approximates the process by using fixed time 
steps, and may lose accuracy as the number of particles grow unless the step size 
is adjusted accordingly. A reasonable way to compare the two methods is thus to 
keep the required accuracy of the tau-leap method fixed. 

From Theorem 4.4 we see that the computational work for the Poisson bridge 
tau-leap method to achieve a relative error e := F,[g(Z T ) — g(Z T )\ , using the time 
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step t = h-y 2 2p , is 

(4.22) ^or/t TL w — w— - 



The comparable work for the SSA is 



WorkssA « — « Ca(Z ) = Ca( 7 Z ) « C 7 P Z P , 



and then we have 



(4.23) « CV"^, 

Thus, asymptotically as 7 — ► 00, for p = 1 the tau-leap method outperforms the 
SSA. For p = 2 the methods are comparable and for p > 2 the SSA seems to be the 
right choice. Note, that for p — 1 the step size r = h is independent of the number 
of particles even though the time between reactions is of order 7~ p . 

Remark 4.13. The estimated work (4.22) and consequently the comparison (4.23) 
is a worst case scenario. For a simple decaying reaction, X — > 0, with propensity 
a(x) = x p for p = 1,2,3, stoichiometric number v = —1 and initial number of 
particles 7, the tau-leap method is 

Xt n+1 =X tn -V(Xf n T n ), p = l,2,3, 
Xt = 7- 

Rewrite 
where 



V(Xf n r n ) = X p Jn + yJXf n T n AW t „ 
V(X p r n )-X p r n 



AW, 



tn ■- 



The Berry-Essen theorem, see [ ], implies that AWj n approaches a N{0,y/f^) dis- 
tributed variable as the number of particles grow. Neglecting the (relatively decreas- 
ing) stochastic term the tau-leap method can thus be approximated by the mean field 
equation 

X' t = -Xf, 4 6(0,71, p= 1,2,3, 
X = 7- 

with scaled solutions 

Zt = e * , Z t = ——r, Z t = 



74+1' ^2724 + 1' 

/orp = 1,2,3, respectively. Since g{Zx) is essentially independent of 7 /orp = 1, 
4/ie worfc is also expected to be independent 0/7, as indicated in (4.22). For p = 2, 3 
the relative solution Zt decays as 7 _1 and most of the decay happen in 7~ p time, so 
for the interval t £ ( 7 ~ p , T) less work is required and using the time-step r — h~f 2 ~ 2p 
in the interval t g (0, 7 ~ p ) implies that the estimate (4.22) gains at least one order 
ofp- 

Remark 4.14. In [2], the authors consider the case of constant density, by assum- 
ing that 

i=l 
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in (2.3). This implies that the propensity can be written as 

a,j(z) = ^(jx) = jA(x), 

for some uniformly bounded function A, and the 'effective ' propensity is thus linear 
in terms ofj. Using this assumption, the authors of[ ] show that choosing t — 7~' 3 , 
for some [3 £ (0, 1), gives a relative error of the same order. Using the same 
constant density assumption would, with our analysis, lead to a step size r = h, 
for some h > 0, and a relative error of the same order, which has the additional 
advantage that the step size is independent 0/7 and the unknown constant j3. 

4.3. Dual approximation. In Theorem 4.4 we showed an a priori estimate for 
the relative global error. For an a posteriori error estimate we have the following 
result: 

Theorem 4.15 (A posteriori error). Assume that Lemma 4-5 and the growth condi- 
tion in Corollary J^.12 hold. Given M u sample paths {X tn (wi)li=l from the Poisson 
bridge tau-leap method, with deterministic time steps {t n }n=o> the relative error can 
be approximated by 



1 

(4.24) E[g(Z T ) - g(Z T )] = — £ /(«,) + O (7"^ ox ) + e d + e u , 

i—l 

where 

N-l , M 

n=0 \j=l 

Here, the time steps are r„ := t n+ \ — t n and the dual weight cpt„ & is defined by 
the backward problem 

(4-25) . + 

<Pt = 9 (Z T ), 

with J tn e R dxd defined by 

(4.26) J tn :=Id+Y^ i>j ■ VS^ZtJ ( r n + — 2==AWi a ] l 5( ^j >0) 
where 

(4-27) AWi n := ^7^^ — , 

are approximate Wiener increments, defined for a{Z trl ) > 0. The process Z tn has 
increments AY^ , as described by (3.6) in Section 3.1. 

Finally, e u is the standard Monte Carlo error with \&r[e u ] = M~ Var[I], and 
Ed is a remainder of diffusion type 

N-l M 
£d = y(^(4+i>Ui) - E[ (p tn+1 I T tn+1 ]j •^2[a J (Zt n+1 ) - a^ZtJ^Vj, 

n=Q 3=1 

where the a -algebra Tt is generated by the history of Z t up to time t. 
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Proof. The proof of the above theorem comes from applying the Trapezoidal rule 
to the error representation formula in Lemma 4.1, Taylor expanding u as in the 
first part of the proof of Theorem 4.4, and using Corollary 4.12. For given time 
steps t„ this gives 



E[g(Z T ) - g(Z T )] 

M N-l , tn+1 

-EE I h 

j=l n=0 
M N-l 



(4.28) 



(aj(Z s ) - dj(Z tn ))(u(Z s + Vj,s) - u(Z s ,s)) 



ds 



= E E T E - aj( z tJ)d z u(Z tn+1 ,t n+1 ) ■ Vj 



j=l n=0 



+ G(7- 1 r„ 3 ) 



M Af-1 



E E f E [(M^ +1 ) -m^jM &„ +1 1 

j=l n=0 



+ 0( 7 - 1 r3)+ ed 



M Af-1 



EE f E fe'(^ +1 )-s-(^j)^„ +1 -^ +0(7-^)+^. 



j=l n=0 



□ 



Remark 4.16. iVoie i/iai (4.26) is the approximation of the first variation dX n+ \/ dX, 
for the chemical Langevin equation, 

M 

(4.29) X tn+1 =X tn +J2 "j (rnCLi (X tn ) + A< yj a, (X tn )) , 

i=l 

where AW 7 ^ ~ iV(0,r„) are standard Wiener increments, and (p tn is thus an ap- 
proximation to the discrete dual for the corresponding value function. This means 
that the Sd term contains errors from both approximating the value function and 
the dual for the Langevin equation. 



5. Examples 

The goal is here to numerically verify the error representation in Lemma 4.1 and 
the a posteriori error estimate in Theorem 4.15. This is a fundamental step for 
future work developing an appropriate adaptive algorithm for the Poisson bridge 
tau-leap method in the spirit of [20]. 

5.1. Testing the estimates. To test the error representation in Lemma 4.1 we 
use the approximation 

N-l M 

E[g(X T ) - g(X T )\ = £ £ E [/,•„] +0 (r^ ax ) , 
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where 

I i,n = ^( a j( x t n +i) -aj{X tn )) 

(5.i) y _ ' _ v 

and show that both Ihs and rhs decay as O (T max ) asymptotically, and that the 
efficiency index Ihs /rhs approaches 1. 

To estimate the true value function u, we solve the Kolmogorov backward equa- 
tion (2.1), which by Assumption 2.1 will become a Ylf = i{Xmax + 1) dimensional 
(stiff) system of ordinary differential equations, since the number of particles in 

(i) 

the chemical system will either be constant or decrease over time, e.g. X t € 
{0,1,..., A^L} for some upper bound X max € Mr. 

For large values of X max the discretization is chosen such that the system is 
solved for a subset of logarithmically distributed integers in [0, Xmax], see Figure 
5.1. Of course, this distribution is in no sense optimal and ideally a spatially 
adaptive algorithm should be used, but it can be expected that points close to any 
XW = will have a greater contribution to the error. 



Logarithmic grid 




Figure 5.1. Example of logarithmic grid used for solving the 
backward Kolmogorov equation in Zjj_. 

Calculating Ihs with sufficient accuracy can be very demanding for problems 
with many species and a high number of particles. A less demanding way is to use 
the approximation 

Ihs = E[g(X T ) - g(X T )} + E[g(X T ) - g(X T )} 
= 2E[g(X T ) -g(X T )} +Q {r 2 max ) , 

"V 

i hs approx 

where the process X t is generated by using half the step size of X t and the sample 
path (A, Y) generated by X t . The estimate lhs approx does not use the value function 
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at all, and has the sample variance of order tM 1 compared to M 1 for Ihs, see 



In practice, the estimate rhs in itself is not much of use since it requires knowl- 
edge of the value function u, and thus rhs must be approximated by computable 
quantities, e.g. as in Theorem 4.15. For this purpose, let 



rhsduai ■= E 



JV-l M 
n=0 j=l 



where the approximate discrete dual is given as in Theorem 4.15 but without the 
scaling factor 7. The accuracy of rhsduai is of great importance to construct a 
proper adaptive algorithm, see e.g. [20, 21, 25]. To show this we use the error 
density 



(5.2) 



Pj,n '■ — 



2r, 



-E[(oj (X n+1 ) - aj (X n )) ip n ■ v 3 



defined for the initial deterministic time steps, i.e. not the time steps given by the 
Poisson bridge tau-leap method in Section 3.1, see [22]. This gives the total error 

JV-l M 

e := rhs dua i = T nPj,n, 

n=0 j=l 

and the current work Work,TL '■— N is then compared with the estimated work to 
achieve the same error e for an optimal adaptive mesh 

/N-l _ \ 2 
Work a = ^ Vf^Tn 



(5.3) 

and a uniform mesh 
(5.4) 



,n=0 



N-l 



Work u = T ^ Pn T n- 



n=0 



5.2. Reactive decay. In this example we have an irreversible reaction where 
molecules of a single species spontaneously disappear, possibly into another particle 
type. The chemical reaction for reactive decay (or isomerization) can be written as 

(5.5) 

and is here described by the linear propensity function 
(5.6) a(x) = cx, 

initial value X(0) = Xq E Z + , and stoichiometric number v = —1. We test four 
different cases: a high number of particles or a low number of particles, that becomes 
negative often or not so often, see Table 5.1. In Figure 5.2 a few realizations of X t 
for t € [0, 1] are shown for each example. 



Example 


1 


2 


3 


4 




10 


10 ti 


10 


10 6 


c 


0.2 


1.0 


2.0 


7.0 



Table 5.1. Reactive decay, cf. (5.5) and (5.6). 
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As quantities of interest we take the first moments of X, i.e. g(x) — x mom 
for mom — 1,2,3 and with final time T = 1. In Figure 5.3, 5.4 and 5.5 the 
convergence of lhs a p prox , rhs and rhsduai can be seen for the different moments. 
In all cases we see a linear 0(t) convergence and in Figure 5.6, 5.7 and 5.8, we see 
that the corresponding efficiency indices stay close to 1, i.e. the computable error 
approximation in Theorem 4.15 agrees well with the error representation formula 
in Lemma 4.1. Comparing the current work W with the work estimates for an 
adaptive mesh in (5.3), and a completely uniform mesh in (5.4), shows that for 
these relatively non-stiff examples adaptivity will not make any improvement, see 
Table 5.2. In this case, a uniform mesh is thus most suitable, and the examples are 
only run for verification purposes. Also note that the error in Example 2 and 4, i.e. 
for a high number of particles, is almost completely governed by the deterministic 
error and require few realizations to achieve a low statistical error of the estimates 
Ihs, rhs and rhsduai- 




Figure 5.2. A few realizations of reactive decay for Example 1 
to 4. These examples are mostly suited for uniform meshes and 
adaptivity will thus not make any improvement. 
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Convergence Convergence 




io _1 10° icr 3 icr 2 io~ 1 

dt dt 



(a) Example 1 (b) Example 2 




Figure 5.3. Reactive decay: Convergence of the errors lhs approx , 
rhs and rhsduai with respect to the maximum time step. Here, 
g(x) = x and the data is scaled by dividing with u(xq, 0). For each 
data point the number of samples is controlled by the standard 
error such that 1.96SE < 0.1/i, where [i is the sample mean and SE 
the standard error of the mean. Bars indicate the 95% confidence 
intervals. 



Example 


1 


2 


3 


4 


Number of realizations 


7.5-10 4 


4.0-10 1 


1.3-10 5 


3.6-10 2 


Current work 


17 


4097 


513 


11265 


Optimal work 


16 


4096 


512 


10866 


Uniform work 


16 


4096 


512 


10868 


Optimal work using dual 


16 


4096 


512 


10867 


Uniform work using dual 


16 


4096 


512 


10869 



Table 5.2. Current and estimated work for first moment. Values 
corresponding to smallest r in Figure 5.3. Two error densities 
are here used: one using the dual as in (5.2) and one using the 



derivative of the true value function. 
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Convergence 



Convergence 



10" 



10" 

































C*dt 










-*-!hs 










. _□_ . r h s 










rhs ^ , 

dual 













10 

dt 

(a) Example 1 
Convergence 




(b) Example 2 
Convergence 




(c) Example 3 (d) Example 4 

FIGURE 5.4. Reactive decay: convergence for second moment of X. 

5.3. Unstable dimer. This stiff model was used in e.g. [1, 13] and has four 
reactions and three species. The reactions are 

Xi->0, X 3 -*2Xi, 

1X\ — > X2, X2 — > X3, 

described by the stoichiometric matrix and propensity function 

'-1 -2 2 \ ( Xl 

1 > = I 1 -1-11 and a(X) = 











1 



0.001JCi(Xi 

0.5X2 
0.04X 2 



1) 



respectively. In Figure 5.9, where a few realizations of the path X t := (X\ , X2 , X3) (t) 
are shown (in a log-lin scale), it can be noted that there is a large difference in 
time-scales; during a very short time most of the Xi-molecules will turn into X2- 
molecules. This difference in time scales has several consequences: the step size 
during the transient phase may have a big effect on the error, and the tau-leap 
method will cause negative populations unless the step size is adjusted accordingly. 

We here let g{x) — X\ + X2 + £3, and the convergence of the corresponding error 
and efficiency index can be seen in Figure 5.10. As in the previous example we see 
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Convergence 



Convergence 




10" 




10" 



10" 
dt 



10 



(a) Example 1 
Convergence 



(b) Example 2 
Convergence 




(c) Example 3 (d) Example 4 

Figure 5.5. Reactive decay: convergence for third moment of X. 

that the error decreases linearly with r, as expected, with error estimates being 
close to each other. 

From Figure 5.11 we see that the error density here play a big role, and that 
the optimal time stepping is to choose very small time steps during the transient 
phase. To see if choosing a pre-leap check as in [12] will give a similar result we 
apply the leap-size condition (3.3) with e = 0.05, which turns out to give almost the 
same error density and proposed time steps. In Table 5.3 a comparison between the 
current and estimated work shows potential for great improvement using adaptive 
time stepping. 



Leap-check 


No 


Yes 


Current work 




2628 


2637 


Optimal work 




573 


575 


Uniform work 




40480 


40501 


Optimal work using 


dual 


529 


531 


Uniform work using 


dual 


42362 


42313 



Table 5.3. Unstable dimer: Current and estimated work for 
realizations. Values corresponding to smallest r in Figure 5.10. 
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Figure 5.6. Reactive decay: Efficiency index for first moment of 
X with 95% confidence interval. The mean and confidence intervals 
were calculated by using 200 Bootstrap samples. 



6. Conclusions 

We have in this work shown a weak global error representation for the tau-leap 
method that can be accurately approximated by a computable leading order term 
using a discrete dual weighted propensity residual, see Lemma 4.1 and Theorem 
4.15. This type of a posteriori error estimates, using discrete dual functions, are 
the first for the tau-leap method and are important tools for the ongoing work on 
developing efficient adaptive time stepping algorithms, see e.g. [22, 21]. Also, we 
have here shown an a priori estimate on the relative error of the tau-leap method, 
that is of order t independently of the number of particles in the system, see 
Theorem 4.4. 

Our results are based on an error representation using the value function for 
the corresponding Kolmogorov backward equation, which for jump processes is 
defined on a discrete lattice. Using extensions from lattices onto real values and 
stochastic representations, we develop weighted estimators for the value function 
and its derivatives, see Lemma 4.5 and Lemma 4.11 respectively. The weighted 
estimators developed here give polynomial bounds on the value function and its 
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Efficiency index 



Efficiency index 




(a) Example 1 
Efficiency index 



(b) Example 2 
Efficiency index 




(c) Example 3 



(d) Example 4 



Figure 5.7. Reactive decay: Efficiency index for second moment 
of X. 



derivatives and improve similar L-infinity bounds in 
the propensity function. 



that are exponential in 
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(a) Convergence 
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Figure 5.10. Unstable dimer: Convergence and efficiency index. 
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Figure 5.11. Unstable dimer: Error densities and step sizes. 
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